#!/usr/bin/python
# Author: Jimmy Saw
# Date: 07-03-2011
# This program parses 454 paired-end reads aligned against an assembled genome.
# Needs MUMMER output file generated like this:
# nucmer genome.fasta 454reads.fasta -p file
# show-coords -r -c -l -T file.delta > file.coords
# Usage: python parse_overlapping_clones.py file.coords > file.matepairs.txt

import sys
import re

mummer_file = sys.argv[1]
mf = open(mummer_file, "rU")
mfl = mf.readlines()

#toprint = int(sys.argv[2])

left = re.compile('(\w+)_L')
right = re.compile('(\w+)_R')

sizemin = 7000
sizemax = 12000

paired = []
unpaired = []
allreads = []
reads_dict = {}
readnames = []

for line in mfl[4:]:
    strand = ""
    l = line.split('\t')
    ref_start = int(l[0])
    ref_stop = int(l[1])
    qry_start = int(l[2])
    qry_stop = int(l[3])
    if qry_stop > qry_start:
        strand = "+"
    else:
        strand = "-"
    per_id = float(l[6])
    qry_cov = float(l[10])
    qry_name = l[12].rstrip()
    readnames.append(qry_name)
    reads_dict[qry_name] = ((ref_start, ref_stop, qry_start, qry_stop, strand, per_id, qry_cov))

mate_pairs = []

for read in readnames:
    if left.match(read):
        pattern = left.match(read)
        mate_pair = pattern.group(1) + "_R"
        if read in reads_dict and mate_pair in reads_dict:
            if reads_dict[read][4] == "+":
                if reads_dict[mate_pair][4] == "-":
                    #calculate distance
                    if reads_dict[mate_pair][1] > reads_dict[read][1]:
                        pair_distance = reads_dict[mate_pair][1] - reads_dict[read][1]
                        ref_start_read = reads_dict[read][0]
                        ref_start_mate = reads_dict[mate_pair][1]
                        mate_pairs.append((ref_start_read, ref_start_mate, pattern.group(1), "L->R", pair_distance))
                    else:
                        pair_distance = reads_dict[read][1] - reads_dict[mate_pair][1]
                        ref_start_mate = reads_dict[mate_pair][0]
                        ref_start_read = reads_dict[read][1]
                        mate_pairs.append((ref_start_mate, ref_start_read, pattern.group(1), "R->L", pair_distance))
                else:
                    #bad alignment. shouldn't be in same orientation as read
                    pass
    elif right.match(read):
        pattern = right.match(read)
        mate_pair = pattern.group(1) + "_L"
        if pattern.group(1) in mate_pairs:
            pass #mate pairs already checked in previous statement
        else:
            if read in reads_dict and mate_pair in reads_dict:
                if reads_dict[read][4] == "+":
                    if reads_dict[mate_pair][4] == "-":
                        #calculate distance
                        if reads_dict[mate_pair][1] > reads_dict[read][1]:
                            pair_distance = reads_dict[mate_pair][1] - reads_dict[read][1]
                            ref_start_read = reads_dict[read][0]
                            ref_start_mate = reads_dict[mate_pair][1]
                            mate_pairs.append((ref_start_read, ref_start_mate, pattern.group(1), "R->L", pair_distance))
                        else:
                            pair_distance = reads_dict[read][1] - reads_dict[mate_pair][1]
                            ref_start_mate = reads_dict[mate_pair][0]
                            ref_start_read = reads_dict[read][1]
                            mate_pairs.append((ref_start_mate, ref_start_read, pattern.group(1), "L->R", pair_distance))
                    else:
                        #bad alignment. shouldn't be in same orientation as read
                        pass
    else:
        #singletons
        #unpaired.append(read)
        pass

mate_pairs.sort()

#Iterate over a given segment of the genome to sort mate pairs
tmpstr = mfl[5].split('\t')
genome_size = int(tmpstr[7])

to_keep_list = []

#Start with mate pairs spanning a given segment of the genome
#Iteration is set to increments of 10,000 bp segments.
i = 0
while i < genome_size:
    read_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+i < start < i+1000 and i+sizemin < stop < i+sizemax:
            read_list.append(pair)
    if len(read_list) > 0:
        #to_keep_list.append(read_list[len(read_list)-1])
        to_keep_list.append(read_list[0]) #only keep one
    i += 10000

#Now, work on pairs that overlap roughly 1000 bp with previous pairs
j = 1000
while j < genome_size - 1000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 2000
while j < genome_size - 2000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 3000
while j < genome_size - 3000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 4000
while j < genome_size - 4000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 5000
while j < genome_size - 5000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 6000
while j < genome_size - 6000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 7000
while j < genome_size - 7000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 8000
while j < genome_size - 8000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 9000
while j < genome_size - 9000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

j = 10000
while j < genome_size - 10000:
    overlap_list = []
    for index, pair in enumerate(mate_pairs):
        start = pair[0]
        stop = pair[1]
        if 0+j < start < j+1000 and j+sizemin < stop < j+sizemax:
            overlap_list.append(pair)
    if len(overlap_list) >= 2:
        #to_keep_list.append(tmp_list[len(tmp_list)-1])
        to_keep_list.append(overlap_list[0]) 
    j += 10000

#Sort the final list according to genome coordinates
to_keep_list.sort()

#Print the mate pairs
for pair in to_keep_list:
    if sizemin < pair[4] < sizemax: #only print pairs with expected insert size
        print pair[2] + "\t" + pair[3] + "\t" + str(pair[0]) + " -> " + str(pair[1]) + "\t" + str(pair[4])

mf.close()
